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Qs . We study the Metropolis dynamics of a simple spin system without dis- 

order, which exhibits glassy dynamics at low temperatures. We use an imple- 
mentation of the algorithm of Bortz, Kalos and Lebowitz jl]. This method 
turns out to be very efficient for the study of glassy systems, which get trapped 

^D ■ in local minima on many different time scales. We find strong evidence of ag- 

. ing effects at low temperatures. We relate these effects to the distribution 

CN , function of the trapping times of single configurations. 
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In the last few years, a lot of work has been devoted to the description of systems with 
glassy dynamics. The phenomenon of aging |^ seems to be ubiquitous. Dynamical aging, 
which may be loosely defined as the property that correlation or response functions on a time 
scale r depend on the age of the system t^ when t c^t^, has been used as a powerful probe in 
experimental studies of the spin glass phase 0. Recently, a lot of numerical studies of various 
spin glass systems have reported aging effects [Q. The theoretical understanding of aging 
is still not satisfactory, although there have been interesting proposals of phenomenological 
models |§|-§1, and also some attempts at a microscopic study and a link with the static 
replica solution of some simple mean field systems ||^r3|. 



In the spin glass case, slow dynamics is generated by the frozen disorder. In this paper 
we want to report on the numerical observation of aging effects in a model which has no 
quenched disorder. This model is a simple spin system with long-range four-spin interactions 
which has been studied from a physical point of view by Bernasconi [^. That quenched 
disorder is not mandatory for the existence of a glassy dynamics is evident: the glassy state 
is a notorious counter example. The system we study provides a simple model for this 
effect. In a sense, the system generates its own disorder through the slowing down of some 
degrees of freedom. This model can help to provide a new bridge between glassy systems 
with or without quenched disorder, on top of the approaches through mode coupling theory 
n|, or some other glass models studied recently |]15||. The emergence of long time scales in 



frustrated systems without disorder has been noticed in recent papers |]16 . 

The numerical observation of some non-equilibrium phenomena which resemble aging in 
a Monte Carlo simulation is rather easy. A more difficult task is to assert whether these 
are just transient phenomena taking place on time scales shorter than some equilibration 
time tf,g or if t^g is infinite. This task is difficult almost by definition: in the glass phase the 
dynamics is slow, and aging should be demonstrated by the existence of relaxation processes 
taking place on very different time scales, in order to be safely distinguished from a simple, 
but slow relaxation rate towards equilibrium. This requires simulating a system at low 
temperatures, in a regime where the acceptance rate of the Metropolis algorithm becomes 
very small. 

This problem of the low acceptance probability can be circumvented efficiently with an 
algorithm, originally proposed by Bortz et ai [|l|, and used in some rare occasions |jl6| [|l^. We 
introduce here an efficient implementation of this algorithm, and apply it to the Bernasconi 
model. At low temperatures, we can simulate this system on time scales which are several 
orders of magnitude larger than those accessible by a straight Monte Carlo simulation. We 
suspect that this algorithm could be quite useful for the simulations of a whole class of glassy 
systems at low temperatures, which we shall characterize briefly in the conclusion. 

The problem studied by Bernasconi some years ago consists in flnding low autocorre- 
lation binary sequences (LABS). This optimization problem, originally introduced for its 
applications in communications, can be turned into a physical problem of N Ising spins, 
where the 2^ spin conflgurations are weighted by a Boltzmann factor, with an energy given 
by: 
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Monte Carlo simulations have been carried out by Bernasconi M, Golay [|18|, and more 



recently by Migliorini [0], who have shown that finding the ground state (or low lying con- 
figurations) is numerically hard; the system freezes at a certain temperature, below which 
the specific heat becomes very small: the system gets trapped into metastable states. Here 
we would like to consider it as a simple prototype of systems exhibiting dynamical aging and 
undergoing a glass transition. Simultaneously to our numerical simulations, several inde- 
pendent works have tried to investigate this problem [^], or another version with periodic 
boundary conditions |2T[], analytically by approximating it by a disordered system. The 
corresponding disordered problem can be studied by the replica technique and undergoes a 
"one step replica symmetry breaking" phase transition characteristic of a class of spin glass 



transitions of the category of the Random Energy Model p2| , |23| . Such a freezing scenario is 
well confirmed by our simulations: at low temperatures, the system gets trapped into single 
configurations which it takes many Monte Carlo steps to escape from. 

Let us now describe the algorithm we have been using. We adopt the usual single spin-fiip 
dynamics, which produces a sequence of configurations 

o"i(t-i) -^ o-2(r2) ... ^ cri{Ti) -^ o-i+i(ri+i) ... etc (2) 

In eq. (|]), the symbols indicate that the system stays in configuration ai for a time Tj etc. 
The transition probability is non-zero only for configurations differing by a single spin. (We 
denote in the following by a^"^^ the configuration which results from a after fiipping spin m). 
We use the Metropolis dynamics, in which the transition probability per unit time At is 
given by the explicit formula: 

^^ ' iV \exp[-/3(E(aH)-E(a))] [otherwise) ^^■■■^'^ l-^J 

Usually, the sequence in eq. (^ is directly simulated by means of a rejection method. 
The move a -^ a^^ is accepted with probability eq. (|]), otherwise the system stays in 
configuration a. The method we use calculates the trapping times Xj not by repeatedly 
rejecting moves, but by sampling Tj from its probability distribution -Po-(t); it then cashes 
in repeatedly on the initial expenditure of calculating Pair). 

When one reaches a configuration a, the transition probabilities p^ = p(c" -^ cr'^J) 
towards a// the neighboring configurations erf™], m = 1,...,N are calculated. The escape 
probability from configuration a is Pesc = J2mPm, and the probability of leaving a after 
exactly r MC is equal to PaiT) = Pesc(l ~ PescY'^"^^ ■ The conditional probability for reaching 
configuration a^"^\ given that one leaves a, is g^ = Pm/Pesc- From a, the algorithm samples 
the time r at which the next visited configuration will be reached, using the distribution 
Po- ('?")• It samples independantly the new configuration a^"^^ (automatically different from 
a) by choosing the value of m with probability g^. This algorithm is exactly equivalent to 
the standard one. 

The price one pays in this algorithm is to compute for each visited configuration all the 
Prni m = 1, ...,N, while only one is computed in the usual method. Therefore this algorithm 
becomes more efficient roughly speaking when the acceptance rate of the usual rejection 
method becomes of order 1/N ^A\. Furthermore we have implemented in the algorithm 



an important feature which provides a large gain in computer time, at least on the LABS 
problem. We constantly keep in memory an "archive" of the different configurations visited 
recently during the simulation and of the corresponding probabilities Pm for m = 1, ..., N. 



Upon encountering a configuration a which belongs to the archive we use the Pm in memory. 
In our implementation the archive is emptied automatically as soon as the number of archived 
configurations reaches a number March- In the actual calculations described below, for the 
runs at the smallest temperature T = .075, a configuration which enters the archive (with 
Narch ~ 400) typically remains active for a very long time (during ~ 10^ configurations), 
before it is moved out. 

Let us mention that the physical limit corresponding to a vanishing elementary time step 
At — i> is trivially taken in this algorithm. In this limit, the distribution of the escape time 
r from the configuration a is 

P,(r) = Aexp[-Ar], (4) 

where A is related to the transition probabilities of (H) by A = Y.rnPi^ ~^ o-[™])/Ar. 

In the above method, the sampling of the dynamical evolution eq. (^ is split into two 
processes, the generation of the trajectories, and the sampling of physical escape times with 
the probability distribution Pa{^)- These two processes are independent sources of statistical 
noise and the second one can be removed. Consider a trajectory (Ti(Ai) -^ cr2{\2) •••—>■ 
cm{^m), where A, is the mean escape time of configuration cxj, characterizing PaiiT) as in 
eq. (^). For any realization, the physical times have to be sampled from eq. (Q). The 
probability Pm{'T') to be at configuration ctm at a given time r is a convolution which can be 
computed by the calculus of residues. In the special case in which all the Aj's are different 
from each other, the result takes the simple form: 

M M \ 

PM{r) = Y. E YtV[1-^^p(-^.^)] (5) 

Unfortunately, we have not succeeded in calculating this quantity in a stable way for large 
M (such as M ~ 10^ ~ 10^), which is the case of practical interest here. A practical answer, 
and the one which we have adopted, is simply to calculate, for a given trajectory eq. (0), a 
fair number of realizations of the escape times (with a completely unsophisticated program, 
we can calculate an average over ~ 100 realizations of the escape times in just about twice 
the time it takes to generate the trajectory). 

We have used this algorithm to study the LABS model. We begin with the results 
from simulated annealing. In fig. 1 we plot the internal energy per spin vs temperature at 
A^ = 100. The annealing scheme starts from a temperature Ti = .6. After thermalizing 
at this temperature, the temperature is lowered down to Tq = .04 using the logarithmic 
scheme: T{t) = Ti — (Ti — TQ)Log{t)/Log{tm), where the duration of the run, tm (measured 
in Monte Carlo steps per Spin (MCS)), takes the values 1.25^0, 1.5^0, 1.752°, 2^°, 2.25^°, 2.5^°. 
So the logarithmic cooling rate dT/d{ln{t)) varies from .125 to .0306. The data is averaged 
over 400 to 50 realizations of the random process of eq. (0) with random initial conditions; 
for each sequence of configurations we have studied a single realization of the escape times. 
The system undergoes a freezing phenomenon at a temperature which depends on the log- 
arithmic cooling rate, reminiscent of the observations in glasses. An extrapolation of the 
zero-temperature energy as a function of the logarithmic cooling rate is compatible with a 
zero-temperature energy per spin of the order of .07 ~ .08. This result is about twice as large 
as the value Ec = .0406. Ec is conjectured by the Bernasconi-Golay approximation ^,TE 



recently derived using an approximate description of this system by one with quenched dis- 



order ^0[ , and is also compatible with exact enumerations of small size systems. We have 
done simulations at A^ = 200, A^ = 400 and N = 401, and t^ = 1.2520,1.520,1.7520,22° 
which are nearly indistinguishable from the data of fig. 1. The difference between the result 
of simulated annealing and Ec may well be due to a dynamic freezing effect, which forbids 
to find the exact ground state for large size systems. 

We have performed a more detailed study of the dynamics at temperatures T = .075, 
which is clearly in the low T regime on all the time scales we can reach, and at T = .25 
for comparison. We start by reviewing the results at the lower temperature, T = .075. 
In fig. 2 we plot the correlation function C{t^ + T.t^) = {l/N)J2i < (^i(tw + T)<^i(tw) > 
vs the time r for various waiting times t^, ranging up to 10^ (from now on all the times 
are measured in MCS). The aging effect is clear, on all the range of waiting times we can 
reach. The behavior is well approximated by a function of r/t^, decreasing as a power law 
(r/t^)~'026 at large arguments, as shown in fig. 3. This is similar to what has been found 
in experiments [0, and numerical simulations of spin glasses 0| or other disordered systems 
p5| . Recent analytical studies of some mean field spin glasses confirm the possibility of such 



a scaling [^ for a category of spin glasses characterized by a first order ( "one step" ) replica 
symmetry breaking scheme like the Random Energy Model P^J23[ ], while the spin glasses 
with full replica symmetry breaking seem to have a more complicated behavior [1^,1^. We 



exhibit here a clear aging effect for the LABS model which is a system without any quenched 
disorder. This data again suggests that the LABS could be a good toy model for the glass 
transition. 

On this aging data one also notices that the asymptotic dynamics is essentially frozen: 
lirriT-^oolif^tni^ooCityj + r, t^) ~ 1. It is interesting to look at a single run of a system with 
N = 400 spins, at T = 0.075. We show such a run, typical of what we observe, in fig. 
4. The simulation, which has taken a few hours on a HP 730 workstation, was stopped 
after a time of r = 10^ MCS. It is immediately seen on the plot of the energy vs time, 
that the energy of the system is not equilibrating on this time scale. Clearly the system 
gets trapped into metastable states with a long lifetime which are well identified by the 
plateaux in the energy. The fluctuations along each plateau are quite small, which is in 
agreement with the above observations: the metastable states at low temperature involve a 
small number of configurations. In contrast we show in figs 5 and 6 the analogous data at a 
higher temperature T = .25. There is no trace of the plateaux in the energy, and the system 
equilibrates after a finite time t ~ 10^, which agrees with the simulated anealing data from 
fig. 1. On the correlation function one observes an interrupted aging effect. 

The trapping behavior is reminiscent of the simple generic model of aging proposed by 
Bouchaud 0. At a qualitative level one does see on fig. 4 that the trapping times (length 
of the plateaux) become longer when the waiting time increases. A quantity which is easily 
accessible with our algorithm is the distribution of escape time from single configurations 
f = 1/A. (Notice that this is not exactly the same as the trapping times because the plateaux 
involve several configurations). In fig. 7 we plot /^Pfc(r') dr' vs r, where Pk{T') is the 
distribution of escape time for the kth visited configuration in the sequence eq. (|). We see 
a clear increase of the typical escape time with k (and therefore with the age of the system). 
For large waiting time (large k) the distribution Pk{'^') is broad, and compatible with the 
idea of M that the average escape time is infinite. The behavior of Pk{T) at large times 



r can be approximated by a power law Pk{T) — t~^'^. The broadness of this distribution 
explains why the algorithm we use is so efficient on this problem. It circumvents the problem 
of trapping by a single configuration, but does not get rid of "renormalized traps" consisting 
of a group of configurations separated from the rest by barriers. A cluster Monte Carlo 



algorithm [^ has been developed to handle this problem. 

To conclude, we have found a model without quenched disorder which is another repre- 
sentative of a well known class of glassy systems, those which undergo a "one step replica 
symmetry breaking" transition in the static approach. Other members of this class include 
the Random Energy Model [^,^, the binary perceptron ||2^, the p-spin spherical spin 



glass [^ , random heteropolymers p9|-pT| , the manifolds in random media with short range 
disorder ||32|, etc. We believe that the algorithm we have described will be efficient for any 



of these problems. Our code is available by e-mail. 
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FIGURE CAPTIONS 

Fig. 1: Simulated annealing results of the energy per spin vs temperature in a iV = 
100 spin system, corresponding to logarithmic cooling rates (from top to bottom) 
.125,.07, .05, .04,.034, .03. 

Fig. 2: The two times correlation of a A^ = 400 spin system at temperature T = .075, 
Citw + T,tyj) vs T, plotted for various values of the waiting time t^ = 5" with n = 
0, 1, ..., 10. The data is averaged over 300 initial conditions (and for each of them 50 
realizations of the sampling of escape times). 

Fig. 3: The same data as fig. 2 plotted vs r/t^, for t^; = 5'^^ with n = 0,1, ..., 5. 

Fig. 4: Energy per spin vs logio(time) in a single run for a A^ = 400 spin system at 
T = .075. 

Fig. 5: The two times correlation of a A^ = 400 spin system at temperature T = .25, 
Citw + T,tw) VS T, plotted for various values of the waiting time t^ = 2.5" with 



n = 0, 1, ..., 7. The data is averaged over 400 initial conditions (and for each of them 
50 reahzations of the samphng of escape times). 

Fig. 6: Energy per spin vs logio(time) in a single run for a iV = 400 spin system at 
T = .25. 

Fig. 7: The integrated distribution of escape times, J^ Pki^') dr', for the kth visited 
configuration, is plotted vs t ioi k = 4" with n = 2, ..., 7, in a simulation of a A^ = 400 
spin system at T = .075. 
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